---
title: "B-MDU_Code"
author: "Kodai Tachibana"
date: "2022/12/2"
output: html_document
editor_options: 
  chunk_output_type: console
---

#package
```{r}
#rm(list=ls())
library(tidyverse)
library(ggplot2)
library(ggrepel)
library(cmdstanr)
library(psych)
library(smacof)
library(bayesplot)
library(rstan)
library(GPArotation)
rstan::rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())
```

#Pre-processing of data loading and labeling
```{r}
mean_ches_data <- read.csv("CHES2019V3.csv")
mean_data <- read.csv("mean_data.csv")
data_ches <- read.csv("CHES2019_experts.csv")
data <- read.csv("raw_data.csv")
#France
#mean_ches_data$party[mean_ches_data$party == "FI"] <- "LFI"
mean_ches_data$party[mean_ches_data$party == "RN"] <- "FN/RN"
mean_ches_data$party[mean_ches_data$party == "LREM"] <- "LaREM"
mean_data <- mean_data %>% 
  dplyr::filter(party != "LP" & party != "NPA" & party != "Rad" & party!= "UDI") 
#data_ches$party_name[data_ches$party_name == "FI"] <- "LFI"
data_ches$party_name[data_ches$party_name == "RN"] <- "FN/RN"
data_ches$party_name[data_ches$party_name == "LREM"] <- "LaREM"
data_ches$party_name[data_ches$party_name == "RPR/UMP"] <- "LR"
data_ches$party_name[data_ches$party_name == "VERTS"] <- "EELV"
data_ches$party_name[data_ches$party_name == "UDF/MODEM"] <- "MoDem"
data <- data %>% dplyr::filter(party != "LP" & party != "NPA" & 
                                 party != "Rad" & party!= "UDI" )

#Germany
mean_ches_data$party[mean_ches_data$party == "GRUNEN"] <- "B90Grune"
mean_data$party[mean_data$party == "Linke"] <- "LINKE"
mean_ches_data <- mean_ches_data %>% dplyr::filter(!(party %in% c( "DieTier", "Piraten")))
data_ches$party_name[data_ches$party_name == "GRUNEN"] <- "B90Grune"
data_ches <- data_ches %>% dplyr::filter(!(party_name %in% c( "DieTier", "Piraten")))

#Italy
mean_ches_data <- mean_ches_data %>% dplyr::filter(party != "SVP" & party != "RI")
mean_ches_data$party[mean_ches_data$party == "LFI"] <- "FI"
mean_data <- mean_data %>% dplyr::filter(party != "MDP")
data_ches <- data_ches %>% dplyr::filter(party_name != "SVP" & party_name != "RI" & party_name != "RAD/R")
mean_ches_data$party[mean_ches_data$party == "LFI"] <- "FI"
data <- data %>% dplyr::filter(party != "MDP")

#Netherland
mean_data$party[mean_data$party == "50Plus"] <- "50PLUS"

```

#Get POPPA solution with smacof to create initial values
```{r}
#-----------------------------------------------------
#Extract country data(France:11,Germany:12,Italy:16,Netherland:19)
nation_data <- mean_data %>%
  dplyr::filter(country_id == 19)

#Extract the columns you need and change the order of the columns(order:lrecon(8),nativism(11),antielitism(5))
nation_select <- nation_data%>%
  dplyr::select(party,lrecon,nativism,antielitism,manichean,indivisble,
                generalwill,peoplecentrism,
                complex,emotional,immigration,eu,
                laworder,lifestyle,intradem,personalised)%>%
  arrange(party)

#---------------------------------------------------------------------------

#invert scale
nation_select$manichean <- 11 - nation_select$manichean
nation_select$indivisble <- 11 - nation_select$indivisble
nation_select$generalwill <- 11 - nation_select$generalwill
nation_select$peoplecentrism <- 11 - nation_select$peoplecentrism
nation_select$antielitism <- 11 - nation_select$antielitism
nation_select$nativism <- 11 - nation_select$nativism
nation_select$laworder <- 11 - nation_select$laworder
nation_select$emotional <- 11 - nation_select$emotional
nation_select$personalised <- 11 - nation_select$personalised
nation_select$intradem <- 11 - nation_select$intradem

#Creating Name Vectors
name_df <- as.data.frame(nation_select$party)
colnames(name_df) <- "party"
#Delete unnecessary column (PARTY)
nation_select_d <- nation_select%>%
  dplyr::select(-party)

#unfolding
set.seed(123)
stressvec_p <- rep(NA,1000)
s1_p <- matrix(NA,nrow = nrow(nation_select_d),ncol = 2)
s2_p <- matrix(NA,nrow = ncol(nation_select_d),ncol = 2)

#Use the least-squares solution as a baseline
fitbest_p <- smacofRect(nation_select_d,
                      ndim = 2,
                      type = "interval")

s1_mean_p = fitbest_p$conf.row 
s2_mean_p = fitbest_p$conf.col

stressvec_p[1] <- fitbest_p$stress
nrun=1000
fit_list_p = vector(mode = "list", length = nrun)
fit_list_p[[1]] = fitbest_p

for(i in 2:nrun){
  s1_p=s1_mean_p + matrix(rnorm(ncol(s1_mean_p)*nrow(s1_mean_p), 0, sd=0.2), 
                      ncol=ncol(s1_mean_p), 
                      nrow=nrow(s1_mean_p))
  s2_p=s2_mean_p + matrix(rnorm(ncol(s2_mean_p)*nrow(s2_mean_p), 0, sd=0.2), 
                      ncol=ncol(s2_mean_p), 
                      nrow=nrow(s2_mean_p))
  fit_list_p[[i]] <- smacofRect(nation_select_d,ndim = 2,type = "interval",
                              init = list(s1_p=s1_p,s2_p=s2_p))
  stressvec_p[i] <- fit_list_p[[i]]$stress
}
#plot(sort(stressvec))
min_num_p = which.min(stressvec_p)
final_init_p = fit_list_p[[min_num_p]]
final_init_p$conf.row
final_init_p$conf.col

#Extracting Data
item_poppa <- final_init_p$conf.col
parties_poppa <- final_init_p$conf.row
item_poppa.dat <- data.frame(
  dim1 = item_poppa[,1],
  dim2 = item_poppa[,2],
  item = as.factor(c(8,11,5,1,2,3,4,6,7,9,10,12,13,14,15))
)
parties_poppa.dat <- data.frame(
  dim1 = parties_poppa[,1],
  dim2 = parties_poppa[,2],
  party = name_df
)
#plot
ggplot()+
  geom_text(data=item_poppa.dat,mapping = aes(x=dim1,y=dim2,label=item),nudge_y = .1,size=4)+
  geom_point(parties_poppa.dat,mapping = aes(x=dim1,y=dim2),size=2.5)+
  geom_text(parties_poppa.dat,mapping = aes(x=dim1,y=dim2,label =party),nudge_y=.1,size=5)+
  theme_bw()+
  labs(x="Dimension 1",y="Dimension 2")+
  theme(aspect.ratio = 1)
```

#Create an observed distance data matrix (POPPA)
```{r}
#Nation:France:11,Germany:12,Italy:16,Netherland:19-----------------
nation_data_indv <- data%>%
  filter(country_id == 19)%>%
  arrange(party)
#---------------------------------------------------------------------------

#Extract the columns
nation_data_individual <- nation_data_indv%>%
  dplyr::select(expert_id,lrecon,nativism,antielitism,manichean,indivisble,
                generalwill,peoplecentrism,
                complex,emotional,immigration,eu,
                laworder,lifestyle,intradem,personalised)

#invert scale
nation_data_individual$manichean <- 11 - nation_data_individual$manichean
nation_data_individual$indivisble <- 11 - nation_data_individual$indivisble
nation_data_individual$generalwill <- 11 - nation_data_individual$generalwill
nation_data_individual$peoplecentrism <- 11 - nation_data_individual$peoplecentrism
nation_data_individual$antielitism <- 11 - nation_data_individual$antielitism
nation_data_individual$nativism <- 11 - nation_data_individual$nativism
nation_data_individual$laworder <- 11 - nation_data_individual$laworder
nation_data_individual$emotional <- 11 - nation_data_individual$emotional
nation_data_individual$personalised <- 11 - nation_data_individual$personalised
nation_data_individual$intradem <- 11 - nation_data_individual$intradem

#party list
party_list_poppa <- nation_data_indv$party%>%
  unique()

#extract expert
id_poppa_list <- nation_data_individual$expert_id%>%
  unique()

#Create a data frame to count the number of item
item_poppa_count <- nation_data_individual%>%
  dplyr::select(-expert_id)

#Create a storage location for "for statemnts"
x_poppa <- array(NA,dim = c(nrow(name_df),#P:政党
                      ncol(item_poppa_count),#J_poppa:項目
                      length(id_poppa_list)))#I_poppa:回答者

#Stores the distance matrix graded by each individual as an array
for(i in 1:length(id_poppa_list)){
  nation_in <- nation_data_individual%>%
    filter(expert_id == id_poppa_list[i])%>%
    dplyr::select(lrecon,nativism,antielitism,manichean,indivisble,
                generalwill,peoplecentrism,
                complex,emotional,immigration,eu,
                laworder,lifestyle,intradem,personalised)%>%
    as.matrix()
  x_poppa[ , , i] <- nation_in
}
#Create an array of 1 for observed, 0 for missing
delta_poppa <- ifelse(is.na(x_poppa),999,x_poppa)
obs_data_poppa <- ifelse(is.na(x_poppa), 0, 1)
```

#bayes set(POPPA)
```{r}
P_poppa <- nrow(name_df) #parties
J_poppa <- ncol(delta_poppa) #items
I_poppa <- nrow(as.data.frame(id_poppa_list)) #experts
K <- 2
unfold_stan_poppa <- list(I=I_poppa,P=P_poppa,K=K,J=J_poppa,delta=delta_poppa,obs_data=obs_data_poppa)
cmdstanpath <- "/Users/tachibananokoudai/.cmdstan/cmdstan-2.27.0"
set_cmdstan_path(cmdstanpath)

model_poppa <- cmdstan_model("B-MDU.stan", quiet=F)
```

#MCMC sampling
```{r}
#Initial value for MCMC
init_poppa <- list(Z_raw=item_poppa,#items
                   X=parties_poppa #parties
                        ) 

unfold_poppa <- model_poppa$sample(data=unfold_stan_poppa,
                               init = function() init_poppa,
                               iter_warmup = 4000,
                               iter_sampling = 10000, 
                               chains = 20,
                               parallel_chains = 20,
                               refresh = 1000,
                               seed = 1)

#Check convergence and log posterior density
unfold_poppa$print("Z",max_rows = 26)
unfold_poppa$draws("lp__") %>% mcmc_trace()
```

#MCMC sampling (only if no convergence, turn again with the largest log posterior density as the initial value)
```{r}
#Examine the chain with the highest posterior log density
jigo <- unfold_poppa$draws("lp__")%>%
  colMeans()%>%
  which.max()

#Take the average of that chain.
X_jigo <- unfold_poppa$draws("X")[,jigo,] %>% apply(3,mean) %>% as.matrix()
Z_raw_jigo <- unfold_poppa$draws("Z_raw")[,jigo,] %>% apply(3,mean) %>% as.matrix()

#Update Initial Values
init_poppa <- list(Z_raw=Z_raw_jigo,
                         X=X_jigo
                         ) 

#sampling
unfold_poppa <- model_poppa$sample(data=unfold_stan_poppa,
                               init = function() init_poppa,
                               iter_warmup = 4000,
                               iter_sampling = 10000, 
                               chains = 4,
                               parallel_chains = 4,
                               adapt_delta = 0.99,
                               max_treedepth = 12,
                               refresh = 1000,
                               seed = 1234
                               )

unfold_poppa$print("Z",max_rows = 26)
unfold_poppa$draws("lp__") %>% mcmc_trace()
```


#target rotation(POPPA)
```{r}
#extract MCMC sample(EAP)
#parties sample
X_poppa <- unfold_poppa$draws("X") %>% apply(3,mean) %>% array(dim=c(P_poppa,K))
#items sample
Z_poppa <- unfold_poppa$draws("Z") %>% apply(3,mean) %>% array(dim=c(J_poppa,K))

#create target matrix for target rotation--------------------------------------
lr_data_poppa <- nation_data%>%
  dplyr::select(party,lroverall)%>%
  arrange(party)%>%
  dplyr::select(-party)

#scaling(Parties)
X_scaling <- cbind((X_poppa[,1]-mean(X_poppa[,1]))/sd(X_poppa[,1]),
                 (X_poppa[,2]-mean(X_poppa[,2]))/sd(X_poppa[,2]))
#Scaling(Items) 
Z_scaling <- cbind((Z_poppa[,1]-mean(X_poppa[,1]))/sd(X_poppa[,1]),
                  (Z_poppa[,2]-mean(X_poppa[,2]))/sd(X_poppa[,2]))

zero_vector_poppa <- as.data.frame(rep(0,nrow(lr_data_poppa)))
target_matrix_poppa <- cbind(lr_data_poppa$lroverall,zero_vector_poppa)%>%
  as.matrix()
#-----------------------------------------------------------------------------------


#Target rotation
X_poppa_targetT <- GPArotation::targetT(L=X_scaling ,Target = target_matrix_poppa) 
X_poppa_rot <- X_poppa_targetT$loadings
#?GPArotation::targetT
#get rotate matrix
rotate_matrix_poppa <- X_poppa_targetT$Th

#Rotate item coordinates by using rotate matrix
Z_poppa_rot <- Z_scaling %*% rotate_matrix_poppa


item_poppa.mc <- data.frame(
  dim1 = Z_poppa_rot[,1],
  dim2 = Z_poppa_rot[,2],
  item = as.factor(c("P8","P11","P5","P1","P2","P3","P4","P6","P7","P9","P10","P12","P13","P14","P15"))
)

parties_poppa.mc <- data.frame(
  dim1 = X_poppa_rot[,1],
  dim2 = X_poppa_rot[,2],
  party = name_df
)

#POPPA_plot
poppa_plot <- ggplot(data = item_poppa.mc,mapping = aes(x=dim1,y=dim2,label=item))+
  ggrepel::geom_text_repel(ggplot2::aes(label = item),color="#B2182B")+
  geom_point(size=2,alpha=0.8,shape = 15,color="#B2182B")+
  geom_text_repel(data = parties_poppa.mc,mapping = aes(x=dim1,y=dim2,
                                                        label = party),size=5)+
  geom_point(data = parties_poppa.mc,mapping = aes(x=dim1,y=dim2,label=party),size=3,alpha=0.8)+
  theme_bw()+
  labs(x="Left-Right",y="Populism")+
  theme(aspect.ratio = 1)+
  labs(title = "POPPA: Netherland")+
  theme(plot.title = element_text(size = 15, hjust = 0.5, face="bold"))+
  labs(x = expression(paste("Left-Right (",italic("r"),"=.83)")),
       y = expression(paste("Populism (", italic("r "), "=.92)")))#+
  #xlim(-2.5,3)+
  #ylim(-2.5,3)
  

#cor
cor_poppa <- nation_data%>%
  arrange(party)%>%
  dplyr::select(lroverall,populism)
cor(X_poppa_rot,cor_poppa)
```



#CHES ver.
```{r}
#-----------------------------------------------------
nation_ches_data <- mean_ches_data %>%
  #dplyr::filter(country == 6)#France
  #dplyr::filter(country == 3)#Germany
  #dplyr::filter(country == 8)#Italy
  dplyr::filter(country == 10)#Netherland

#Eliminate entries containing unnecessary columns (country,salience,dissent,require)
nation_ches_ex <- nation_ches_data%>%
  dplyr::select(!ends_with(c("country","eastwest","salience","dissent","require")))%>%
  arrange(party)

#---------------------------------------------------------------------------
#extract columns and fix order(regions(15),deregulation(6),sociallifestyile(9))
nation_ches_select <- nation_ches_ex%>%
  dplyr::select(party,immigrate_policy:regions)%>%
  relocate(regions, deregulation, sociallifestyle, .after = party)

#invert scale
nation_ches_select$immigrate_policy <- 11 - nation_ches_select$immigrate_policy
nation_ches_select$multiculturalism <- 11 - nation_ches_select$multiculturalism
nation_ches_select$redistribution <- 11 - nation_ches_select$redistribution
nation_ches_select$environment <- 11 - nation_ches_select$environment
nation_ches_select$spendvtax <- 11 - nation_ches_select$spendvtax
nation_ches_select$deregulation <- 11 - nation_ches_select$deregulation
nation_ches_select$econ_interven <- 11 - nation_ches_select$econ_interven
nation_ches_select$civlib_laworder <- 11 - nation_ches_select$civlib_laworder
nation_ches_select$sociallifestyle <- 11 - nation_ches_select$sociallifestyle
nation_ches_select$religious_principles <- 11 - nation_ches_select$religious_principles
nation_ches_select$ethnic_minorities <- 11 - nation_ches_select$ethnic_minorities
nation_ches_select$nationalism <- 11 - nation_ches_select$nationalism
nation_ches_select$urban_rural <- 11 - nation_ches_select$urban_rural
nation_ches_select$protectionism <- 11 - nation_ches_select$protectionism


name_ches_df <- as.data.frame(nation_ches_select$party)
#name_ches_df$`nation_ches_select$party`[name_ches_df$`nation_ches_select$party` == "FI"] <- "LFI"
#name_ches_df$party[name_ches_df$party == "LFI"] <- "FI"
colnames(name_ches_df) <- "party"

nation_ches_select_d <- nation_ches_select%>%
  dplyr::select(-party)
#unfolding
set.seed(123)
stressvec <- rep(NA,1000)
s1 <- matrix(NA,nrow = nrow(nation_ches_select_d),ncol = 2)
s2 <- matrix(NA,nrow = ncol(nation_ches_select_d),ncol = 2)

#baseline
fitbest <- smacofRect(nation_ches_select_d,
                      ndim = 2,
                      type = "interval")

s1_mean = fitbest$conf.row 
s2_mean = fitbest$conf.col

stressvec[1] <- fitbest$stress
nrun=1000
fit_list = vector(mode = "list", length = nrun)
fit_list[[1]] = fitbest

for(i in 2:nrun){
  s1=s1_mean + matrix(rnorm(ncol(s1_mean)*nrow(s1_mean), 0, sd=0.2), 
                      ncol=ncol(s1_mean), 
                      nrow=nrow(s1_mean))
  s2=s2_mean + matrix(rnorm(ncol(s2_mean)*nrow(s2_mean), 0, sd=0.2), 
                      ncol=ncol(s2_mean), 
                      nrow=nrow(s2_mean))
  fit_list[[i]] <- smacofRect(nation_ches_select_d,ndim = 2,type = "interval",
                              init = list(s1=s1,s2=s2))
  stressvec[i] <- fit_list[[i]]$stress
}
#plot(sort(stressvec))
min_num = which.min(stressvec)
final_init = fit_list[[min_num]]
final_init$conf.row
final_init$conf.col

#extract
item_ches <- final_init$conf.col
parties_ches <- final_init$conf.row
item.ches_dat <- data.frame(
  dim1 = item_ches[,1],
  dim2 = item_ches[,2],
  item = as.factor(c(15,6,9,1,2,3,4,5,7,8,10,11,12,13,14))
)


parties.ches_dat <- data.frame(
  dim1 = parties_ches[,1],
  dim2 = parties_ches[,2],
  party = name_ches_df
)
#plot
ggplot()+
  geom_text(data=item.ches_dat,mapping = aes(x=dim1,y=dim2,label=item),nudge_y = .1,size=4,color="gray")+
  geom_point(parties.ches_dat,mapping = aes(x=dim1,y=dim2),size=2.5)+
  geom_text(parties.ches_dat,mapping = aes(x=dim1,y=dim2,label =party),nudge_y=.1,size=5)+
  theme_bw()+
  labs(x="Dimension 1",y="Dimension 2")+
  theme(aspect.ratio = 1)
```

#CHES ver
```{r}
#-----------------------------------------------------
nation_data_ches_indv <- data_ches%>%
  #filter(party_id >= 601 & party_id <= 628) #France
  #filter(party_id >= 301 & party_id <= 312)#Germany
  #filter(party_id >= 811 & party_id <= 845)#Italy
  filter(party_id >= 1001 & party_id <= 1051) #Netherlands
#---------------------------------------------------------------------------

nation_ches_ex_indv <- nation_data_ches_indv%>%
  dplyr::select(!ends_with(c("country","eastwest","salience","dissent","require")))%>%
  arrange(party_name)
nation_data_ches_individual <- nation_ches_ex_indv%>%
  dplyr::select(id,party_name,immigrate_policy:regions)%>%
  relocate(regions, deregulation, sociallifestyle, .after = party_name)


nation_data_ches_individual$immigrate_policy <- 11 - nation_data_ches_individual$immigrate_policy
nation_data_ches_individual$multiculturalism <- 11 - nation_data_ches_individual$multiculturalism
nation_data_ches_individual$redistribution <- 11 - nation_data_ches_individual$redistribution
nation_data_ches_individual$environment <- 11 - nation_data_ches_individual$environment
nation_data_ches_individual$spendvtax <- 11 - nation_data_ches_individual$spendvtax
nation_data_ches_individual$deregulation <- 11 - nation_data_ches_individual$deregulation
nation_data_ches_individual$econ_interven <- 11 - nation_data_ches_individual$econ_interven
nation_data_ches_individual$civlib_laworder <- 11 - nation_data_ches_individual$civlib_laworder
nation_data_ches_individual$sociallifestyle <- 11 - nation_data_ches_individual$sociallifestyle
nation_data_ches_individual$religious_principles <- 11 - nation_data_ches_individual$religious_principles
nation_data_ches_individual$ethnic_minorities <- 11 - nation_data_ches_individual$ethnic_minorities
nation_data_ches_individual$nationalism <- 11 - nation_data_ches_individual$nationalism
nation_data_ches_individual$urban_rural <- 11 - nation_data_ches_individual$urban_rural
nation_data_ches_individual$protectionism <- 11 - nation_data_ches_individual$protectionism

party_list <- nation_data_ches_individual$party_name%>%
  unique()

id_ches_list <- nation_data_ches_individual$id%>%
  unique()

item_ches_count <- nation_data_ches_individual%>%
  dplyr::select(-c(id,party_name))

x_ches <- array(NA,dim = c(nrow(name_ches_df),#P:party
                      ncol(item_ches_count),#J:item
                      length(id_ches_list)))#I:expert


for(i in 1:length(id_ches_list)){
  nation_in <- nation_data_ches_individual%>%
    dplyr::filter(id == id_ches_list[i])%>%
    dplyr::select(-c(id,party_name))%>%
    as.matrix()
  x_ches[ , , i] <- nation_in 
}

delta_ches <- ifelse(is.na(x_ches),999,x_ches)
obs_data_ches <- ifelse(is.na(x_ches), 0, 1)
```

#bayes set(CHES)
```{r}
P_ches <- nrow(name_ches_df) #parties
J_ches <- ncol(item_ches_count) #items
I_ches <- nrow(as.data.frame(id_ches_list)) #experts
K <- 2
unfold_stan_ches <- list(I=I_ches,P=P_ches,K=K,J=J_ches,delta=delta_ches,obs_data=obs_data_ches)
cmdstanpath <- "/Users/tachibananokoudai/.cmdstan/cmdstan-2.27.0"
set_cmdstan_path(cmdstanpath)

model_ches <- cmdstan_model("B-MDU.stan", quiet=F)
```

#MCMC sampling(CHES)
```{r}
#Initial value for MCMC
init_ches <- list(Z_raw=item_ches,
                         X=parties_ches  
                         ) 
#samoling
unfold_ches <- model_ches$sample(data=unfold_stan_ches,
                               init = function() init_ches,
                               iter_warmup = 4000,
                               iter_sampling = 10000, 
                               chains = 30,
                               parallel_chains = 30,
                               adapt_delta = 0.99,
                               max_treedepth = 12,
                               refresh = 1000,
                               seed = 1234
                               )

#Check convergence and log posterior density
unfold_ches$print("Z",max_rows = 26)
unfold_ches$draws("lp__") %>% mcmc_trace()
```


```{r}
#Examine the chain with the highest posterior log density
jigo <- unfold_ches$draws("lp__")%>%
  colMeans()%>%
  which.max()

#Take the average of that chain.
X_jigo <- unfold_ches$draws("X")[,jigo,] %>% apply(3,mean) %>% as.matrix()
Z_raw_jigo <- unfold_ches$draws("Z_raw")[,jigo,] %>% apply(3,mean) %>% as.matrix()

#Update Initial Values
init_ches <- list(Z_raw=Z_raw_jigo,
                         X=X_jigo
                         ) 

#sampling
unfold_ches <- model_ches$sample(data=unfold_stan_ches,
                               init = function() init_ches,
                               iter_warmup = 4000,
                               iter_sampling = 10000, 
                               chains = 4,
                               parallel_chains = 4,
                               adapt_delta = 0.99,
                               max_treedepth = 12,
                               refresh = 1000,
                               seed = 1234
                               )

#Check convergence and log posterior density
unfold_ches$print("Z",max_rows = 26)
unfold_ches$draws("lp__") %>% mcmc_trace()
unfold_ches$save_output_files("cmdstanr")
```

#target rotation(CHES)
```{r}
#extract MCMC samples(EAP)
#parties
X_ches <- unfold_ches$draws("X") %>% apply(3,mean) %>% array(dim=c(P_ches,K))
#items
Z_ches <- unfold_ches$draws("Z") %>% apply(3,mean) %>% array(dim=c(J_ches,K))

#target rotation--------------------------------------
lr_data_ches <- nation_ches_data%>%
  dplyr::select(party,lrgen)%>%
  arrange(party)%>%
  dplyr::select(-party)

#scaling(Parties)
X_scaling <- cbind((X_ches[,1]-mean(X_ches[,1]))/sd(X_ches[,1]),
                 (X_ches[,2]-mean(X_ches[,2]))/sd(X_ches[,2]))
#Scaling(Items) 
Z_scaling <- cbind((Z_ches[,1]-mean(X_ches[,1]))/sd(X_ches[,1]),
                  (Z_ches[,2]-mean(X_ches[,2]))/sd(X_ches[,2]))

#0
zero_vector_ches <- as.data.frame(rep(0,nrow(lr_data_ches)))
target_matrix_ches <- cbind(lr_data_ches$lrgen,
                            zero_vector_ches)%>%
  as.matrix()
#-----------------------------------------------------------------------------------
#Target rotation
X_ches_targetT <- GPArotation::targetT(L=X_scaling ,Target = target_matrix_ches) 
X_ches_rot <- X_ches_targetT$loadings
#Obtaining a rotation matrix
rotate_matrix_ches <- X_ches_targetT$Th

#The rotation matrix is used to rotate the item coordinates by the same amount.
Z_ches_rot <- Z_scaling %*% rotate_matrix_ches

item_ches.mc<- data.frame(
  dim1 = Z_ches_rot[,1],
  dim2 = Z_ches_rot[,2],
  item = as.factor(c("C15","C6","C9","C1","C2","C3","C4","C5","C7","C8","C10","C11","C12","C13","C14"))
)



parties_ches.mc <- data.frame(
  dim1 = X_ches_rot[,1],
  dim2 = X_ches_rot[,2],
  party = name_ches_df
)

#CHES plot
ches_plot <- ggplot(data = item_ches.mc,mapping = aes(x=dim1,y=dim2,label=item))+
  ggrepel::geom_text_repel(ggplot2::aes(label = item),color="#2166AC")+
  geom_point(size=2,alpha=0.8,shape = 4,color="#2166AC")+
  geom_text_repel(data = parties_ches.mc,mapping = aes(x=dim1,y=dim2,label = party)
            ,size=5)+
  geom_point(data = parties_ches.mc,mapping = aes(x=dim1,y=dim2,label=party),size=3,alpha=0.8)+
  theme_bw()+
  labs(x="Left-Right",y="Populism")+
  theme(aspect.ratio = 1)+
  labs(title = "CHES: Netherland")+
  theme(plot.title = element_text(size = 15, hjust = 0.5, face="bold"))+
  labs(x = expression(paste("Left-Right (",italic("r"),"=.63)")),
       y = expression(paste("Populism (", italic("r "), "=.51)")))+
  xlim(-2,3)+
  ylim(-2.5,2.5)
ches_plot
  


#cor
cor_ches <- nation_ches_data%>%
  arrange(party)%>%
  dplyr::select(lrgen,people_vs_elite,members_vs_leadership)
cor(X_ches_rot,cor_ches)
gridExtra::grid.arrange(poppa_plot, ches_plot, ncol = 2)
```


#bayes set(Integrated Model)
```{r}
P <- nrow(name_ches_df) #It does not matter which party (name_df or name_ches_df) is used.
J_ches <- ncol(delta_ches) #items
J_poppa <- ncol(delta_poppa)
I_ches <- nrow(as.data.frame(id_ches_list)) #experts(CHES)
I_poppa <- nrow(as.data.frame(id_poppa_list))#experts(POPPA)
K <- 2
unfold_stan_integration <- list(I=I_ches,I_p=I_poppa,
                                J=J_ches,J_p=J_poppa,
                                P=P,
                                K=K,
                                delta=delta_ches,delta_p=delta_poppa,
                                obs_data=obs_data_ches,
                                obs_data_p=obs_data_poppa)
cmdstanpath <- "/Users/tachibananokoudai/.cmdstan/cmdstan-2.27.0"
set_cmdstan_path(cmdstanpath)

model_integration <- cmdstan_model("Integration_Model.stan", quiet=F)
```

#MCMC sampling(Integration Model)
```{r}
#initial value
init.integration <- list(X=scale(parties_poppa))

#sampling(30chains)
unfold_integration <- model_integration$sample(data=unfold_stan_integration,
                               init = function() init.integration,
                               iter_warmup = 4000,
                               iter_sampling = 10000, 
                               chains = 30,
                               parallel_chains = 30,
                               adapt_delta = 0.99,
                               max_treedepth = 15,
                               refresh = 1000,
                               seed = 1234
                               )

#Check convergence and log posterior density
unfold_integration$print("Z",max_rows = 26)
unfold_integration$draws("lp__") %>% mcmc_trace()
```


#MCMC sampling (only if no convergence, turn again with the largest log posterior density as the initial value)
```{r}
#Examine the chain with the highest posterior log density
jigo <- unfold_integration$draws("lp__")%>%
  colMeans()%>%
  which.max()

#Take the average of that chain.
X_jigo <- unfold_integration$draws("X")[,jigo,] %>% apply(3,mean) %>% as.matrix()
Z_raw_jigo <- unfold_integration$draws("Z_raw")[,jigo,] %>% apply(3,mean) %>% as.matrix()
Z_raw_p_jigo <- unfold_integration$draws("Z_raw_p")[,jigo,] %>% apply(3,mean) %>% as.matrix()

#Update Initial Values
init.integration <- list(X=X_jigo,
                         Z_raw=Z_raw_jigo,
                         Z_raw_p=Z_raw_p_jigo
                         )


#4chains
unfold_integration <- model_integration$sample(data=unfold_stan_integration,
                               init = function() init.integration,
                               iter_warmup = 4000,
                               iter_sampling = 100000, 
                               chains = 4,
                               parallel_chains = 4,
                               adapt_delta = 0.99,
                               max_treedepth = 15,
                               thin = 100,
                               refresh = 1000,
                               seed=1234
                               )

#Check convergence and log posterior density
unfold_integration$save_output_files("cmdstanr")
unfold_integration$print("Z",max_rows = 26)
unfold_integration$draws("lp__") %>% mcmc_trace()
unfold_integration$draws("X[1,1]") %>% mcmc_trace()
```

#target rotation (Integration Model)
```{r}
#Extract MCMC sample(EAP)
#Parties sample
X_integration <- unfold_integration$draws("X") %>% apply(3,mean) %>% array(dim=c(K,P))%>%
  t()
#Items sample
Z_integration_ches <- unfold_integration$draws("Z") %>% apply(3,mean) %>% array(dim=c(K,J_ches))%>%
  t()
Z_integration_poppa <- unfold_integration$draws("Z_p") %>% apply(3,mean) %>% array(dim=c(K,J_poppa))%>%
  t()


#Creating target matrix for Target rotation--------------------------------------
lr_data_integration <- nation_data%>%
  dplyr::select(party,lroverall)%>%
  arrange(party)%>%
  dplyr::select(-party)

#scaling(Parties)
X_scaling <- cbind((X_integration[,1]-mean(X_integration[,1]))/sd(X_integration[,1]),
                 (X_integration[,2]-mean(X_integration[,2]))/sd(X_integration[,2]))
#Scaling(Items) 
Z_scaling_ches <- cbind((Z_integration_ches[,1]-mean(X_integration[,1]))/sd(X_integration[,1]),
                 (Z_integration_ches[,2]-mean(X_integration[,2]))/sd(X_integration[,2]))

Z_scaling_poppa <- cbind((Z_integration_poppa[,1]-mean(X_integration[,1]))/sd(X_integration[,1]),
                  (Z_integration_poppa[,2]-mean(X_integration[,2]))/sd(X_integration[,2]))

#target matrix
zero_vector_integration <- as.data.frame(rep(0,nrow(lr_data_integration)))
target_matrix_integration <- cbind(lr_data_integration,zero_vector_integration)%>%
  as.matrix()
#-----------------------------------------------------------------------------------
#Target rotation
X_integration_targetT <- GPArotation::targetT(L=X_scaling ,Target = target_matrix_integration) 
X_integration_rot <- X_integration_targetT$loadings
#obtain rotate matrix
rotate_matrix_integration <- X_integration_targetT$Th

#rotate item coordinate by using rotate matrix
Z_integration_rot_ches <- Z_scaling_ches %*% rotate_matrix_integration
Z_integration_rot_poppa <- Z_scaling_poppa %*% rotate_matrix_integration

#data frame(Items of CHES)
item.ches_integration <- data.frame(
  dim1 = Z_integration_rot_ches[,1],
  dim2 = Z_integration_rot_poppa[,2],
  item = as.factor(c("C15","C6","C9","C1","C2","C3","C4","C5","C7","C8","C10","C11","C12","C13","C14")))

#data frame(Items of POPPA)
item.poppa_integration <- data.frame(
  dim1 = Z_integration_rot_poppa[,1],
  dim2 = Z_integration_rot_poppa[,2],
  item = as.factor(c("P8","P11","P5","P1","P2","P3","P4","P6","P7","P9","P10","P12","P13","P14","P15")))


#data frame(Parties)
parties.mc_integration <- data.frame(
  dim1 = X_integration_rot[,1],
  dim2 = X_integration_rot[,2],
  party = name_ches_df 
)

#Italy_ches <- item.ches_integration
#Italy_poppa <- item.poppa_integration
#Italy_parties <- parties.mc_integration
#France_ches <- item.ches_integration
#France_poppa <- item.poppa_integration
#France_parties <- parties.mc_integration
#Germany_ches <- item.ches_integration
#Germany_poppa <- item.poppa_integration
#Germany_parties <- parties.mc_integration
#Netherland_ches <- item.ches_integration
#Netherland_poppa <- item.poppa_integration
#Netherland_parties <- parties.mc_integration



library(ggtext)
#plot
ggplot()+
  geom_text_repel(data=item.ches_integration,mapping=aes(x=dim1,y=dim2,label=item),
                  size=4,color="#2166AC",nudge_y = -.1)+
  geom_text_repel(data=item.poppa_integration,mapping=aes(x=dim1,y=dim2,label=item),
                  size=4,color="#B2182B")+
  geom_text_repel(data=parties.mc_integration,mapping=aes(x=dim1,y=dim2,label=party),
                  size=5)+
  geom_point(item.ches_integration,mapping = aes(x=dim1,y=dim2),size=2,alpha=0.8,color="#2166AC",
             shape = 4)+
  geom_point(item.poppa_integration,mapping = aes(x=dim1,y=dim2),size=2,alpha=0.8,color="#B2182B",
             shape = 15)+
  geom_point(parties.mc_integration,mapping = aes(x=dim1,y=dim2),size=3,alpha=0.8)+
  theme_bw()+
  labs(x = expression(paste("Left-Right (",italic("r"),"=.77)")),
       y = expression(paste("Populism (", italic("r "), "=.59)")))+
  theme(aspect.ratio = 1)+
  labs(title = "Integrated model (POPPA & CHES): France")+
  theme(plot.title = element_text(size = 15, hjust = 0.5, face="bold"))


#Correlation
#CHES
cor_df <- nation_ches_data%>%
  dplyr::select(party,lrgen,galtan,people_vs_elite,members_vs_leadership)%>%
  arrange(party)%>%
  dplyr::select(-party)%>%
  scale()
#POPPA
cor_df_p <- nation_data%>%
  dplyr::select(party,lroverall,populism)%>%
  arrange(party)%>%
  dplyr::select(-party)%>%
  scale()
#bind data
cor_data <- cbind(cor_df,cor_df_p)

cor(X_integration_rot,cor_data)
d <- Germany_parties%>% select(-party) %>% as.matrix()
cor(d,cor_data)
```
















